Conducting Meta-Analytic Structural Equation Modeling with R

Mike Cheung

2024 ISDSA Meeting online workshop (21 July 2024)

Introduction

Who is the instructor?

Materials

## Required packages in this workshop
lib2install <- c("metaSEM", "lavaan", "semPlot", "readxl")

## Install them automatically if they are not available on your computer
for (i in lib2install) {
  if (!(i %in% rownames(installed.packages()))) install.packages(i)
}

Objectives of the workshop

A crash course in SEM

What is SEM?

Key concept 1: Observed vs. latent variables

An example

Key concept 2: Models of the relationship among the constructs

Key concept 3: Notations for path diagrams

Path analysis

Confirmatory factor analysis (CFA)

SEM

How to specify a structural equation model?

Model specification

Model-identification (1)

Model-identification (2)

Model-identification (3)

Data in the analysis

Sample data
x1 y
7.49 5.27
10.66 5.29
9.61 4.40
14.43 7.22
10.58 5.37
11.59 4.91
7.09 3.83
13.57 6.53
5.87 2.26
8.20 11.93

Model specification

Parameter estimation (1)

Parameter estimation (2)

Model testing

Model evaluation in SEM

Chi-square test statistics

Issues of chi-square test statistics

Goodness-of-fit indices

Incremental fit indices

Residual-based indices (1)

Residual-based indices (2)

What do we need to report in research articles?

A crash course in meta-analysis

Are systematic review and meta-analysis the same?

Not all types of evidence are equal.

Problems in empirical research

Problems of individual study

Objectives of a meta-analysis

Effect sizes and their sampling variances

What is an effect size? (1)

What is an effect size? (2)

Common families of effect sizes

Cheung and Vijayakumar (2016, Table 1)
Cheung and Vijayakumar (2016, Table 1)

Correlation (r)

Models for meta-analysis

Random-effects model

Differences between a fixed-effect and random-effects model

Fixed-effect model

Borenstein et al. (2009, Figure 11.1)
Borenstein et al. (2009, Figure 11.1)
Borenstein et al. (2009, Figure 11.3)
Borenstein et al. (2009, Figure 11.3)

Random-effects model

Borenstein et al. (2009, Figure 12.1)
Borenstein et al. (2009, Figure 12.1)
Borenstein et al. (2009, Figure 12.2)
Borenstein et al. (2009, Figure 12.2)
Borenstein et al. (2009, Figure 12.4)
Borenstein et al. (2009, Figure 12.4)

Issues in interpreting \(\tau^2\) (1)

Issues in interpreting \(\tau^2\) (2)

Confidence interval

Why are researchers concerned about the heterogeneity of studies?

Exploring heterogeneity

Subgroup analysis

Mixed-effects model

What is MASEM?

MASEM

Benefits of MASEM to SEM

Benefits of MASEM to meta-analysis

Problems in primary research using SEM

Terms used in the literature

Example: Brown and Stayman (1992)

Brown and Stayman (1992)
Brown and Stayman (1992)

Example: Premack and Hunter (1988)

Premack and Hunter (1988)
Premack and Hunter (1988)

Example: Norton et al. (2013)

Norton et al. (2013)
Norton et al. (2013)

Approaches to MASEM

Jak and Cheung (2020, Table 1)
Jak and Cheung (2020, Table 1)

Univariate approach (Viswesvaran & Ones, 1995)

Problems of the univariate approach

One (sample) size does not fit all!

Problems of treating a correlation matrix as a covariance matrix

A simulation study comparing some of these methods

Jak and Cheung (2020)
Jak and Cheung (2020)
Jak and Cheung (2020)
Jak and Cheung (2020)

TSSEM approach

Random-effects TSSEM

A conceptual overview of TSSEM

Jak and Cheung (2020, Figure 1)
Jak and Cheung (2020, Figure 1)

Random-effects TSSEM: Stage 1 analysis (1)

Random-effects TSSEM: Stage 1 analysis (2)

Random-effects TSSEM: Stage 2 analysis

OSMASEM approach

Jak and Cheung (2020)
Jak and Cheung (2020)

Examples and applications

An example: A Five-factor model

The dataset

Reading the Excel file into R

## Load the library for MASEM
library(metaSEM)

## Load the library to read XLSX file
library(readxl)

## Read the study characteristics
study <- read_xlsx("Digman97.xlsx", sheet="Info")

## Display a few studies
head(study)
## # A tibble: 6 × 3
##   Study                               n Cluster     
##   <chr>                           <dbl> <chr>       
## 1 Digman 1 (1994)                   102 Children    
## 2 Digman 2 (1994)                   149 Children    
## 3 Digman 3 (1963c)                  334 Children    
## 4 Digman & Takemoto-Chock (1981b)   162 Children    
## 5 Graziano & Ward (1992)             91 Adolescents 
## 6 Yik & Bond (1993)                 656 Young adults
## Create an empty list to store the correlation matrices
Digman97.data <- list()
  
## Read 1 to 14 correlation matrices
for (i in 1:14) {
  ## Read each sheet and convert it into a matrix
  mat <- as.matrix(read_xlsx("Digman97.xlsx", sheet=paste0("Study ", i)))
  ## Add the row names
  rownames(mat) <- colnames(mat)
  ## Save it into a list
  Digman97.data[[i]] <- mat
}

## Add the names of the studies
names(Digman97.data) <- study$Study

## Show the first few studies
head(Digman97.data)
## $`Digman 1 (1994)`
##        A     C   ES     E    I
## A   1.00  0.62 0.41 -0.48 0.00
## C   0.62  1.00 0.59 -0.10 0.35
## ES  0.41  0.59 1.00  0.27 0.41
## E  -0.48 -0.10 0.27  1.00 0.37
## I   0.00  0.35 0.41  0.37 1.00
## 
## $`Digman 2 (1994)`
##        A    C   ES     E     I
## A   1.00 0.39 0.53 -0.30 -0.05
## C   0.39 1.00 0.59  0.07  0.44
## ES  0.53 0.59 1.00  0.09  0.22
## E  -0.30 0.07 0.09  1.00  0.45
## I  -0.05 0.44 0.22  0.45  1.00
## 
## $`Digman 3 (1963c)`
##       A     C   ES     E    I
## A  1.00  0.65 0.35  0.25 0.14
## C  0.65  1.00 0.37 -0.10 0.33
## ES 0.35  0.37 1.00  0.24 0.41
## E  0.25 -0.10 0.24  1.00 0.41
## I  0.14  0.33 0.41  0.41 1.00
## 
## $`Digman & Takemoto-Chock (1981b)`
##        A     C   ES     E     I
## A   1.00  0.65 0.70 -0.26 -0.03
## C   0.65  1.00 0.71 -0.16  0.24
## ES  0.70  0.71 1.00  0.01  0.11
## E  -0.26 -0.16 0.01  1.00  0.66
## I  -0.03  0.24 0.11  0.66  1.00
## 
## $`Graziano & Ward (1992)`
##       A    C   ES    E    I
## A  1.00 0.64 0.35 0.29 0.22
## C  0.64 1.00 0.27 0.16 0.22
## ES 0.35 0.27 1.00 0.32 0.36
## E  0.29 0.16 0.32 1.00 0.53
## I  0.22 0.22 0.36 0.53 1.00
## 
## $`Yik & Bond (1993)`
##       A    C   ES    E    I
## A  1.00 0.66 0.57 0.35 0.38
## C  0.66 1.00 0.45 0.20 0.31
## ES 0.57 0.45 1.00 0.49 0.31
## E  0.35 0.20 0.49 1.00 0.59
## I  0.38 0.31 0.31 0.59 1.00
## Extract the sample sizes
Digman97.n <- study$n
Digman97.n
##  [1]  102  149  334  162   91  656   70   70  277  227 1000  227   91 1040
## Extract the cluster
Digman97.cluster <- study$Cluster
Digman97.cluster
##  [1] "Children"      "Children"      "Children"      "Children"     
##  [5] "Adolescents"   "Young adults"  "Young adults"  "Young adults" 
##  [9] "Mature adults" "Mature adults" "Mature adults" "Mature adults"
## [13] "Mature adults" "Mature adults"
## Display the no. of studies
pattern.na(Digman97.data, show.na=FALSE)
##     A  C ES  E  I
## A  14 14 14 14 14
## C  14 14 14 14 14
## ES 14 14 14 14 14
## E  14 14 14 14 14
## I  14 14 14 14 14
## Display the cumulative sample sizes
pattern.n(Digman97.data, Digman97.n)
##       A    C   ES    E    I
## A  4496 4496 4496 4496 4496
## C  4496 4496 4496 4496 4496
## ES 4496 4496 4496 4496 4496
## E  4496 4496 4496 4496 4496
## I  4496 4496 4496 4496 4496

Proposed model

model <- "## Factor loadings
          ## Alpha is measured by A, C, and ES
          Alpha =~ A + C + ES
          ## Beta is measured by E and I
          Beta =~ E + I
          ## Factor correlation between Alpha and Beta
          Alpha ~~ Beta"

## Display the model
plot(model, color="yellow")

## Convert the lavaan syntax into a RAM model as the metaSEM only knows the RAM model
## It is important to ensure that the variables are arranged in A, C, ES, E, and I.
RAM <- lavaan2RAM(model, obs.variables=c("A","C","ES","E","I"))
RAM
## $A
##       A   C   ES  E   I   Alpha           Beta         
## A     "0" "0" "0" "0" "0" "0.1*AONAlpha"  "0"          
## C     "0" "0" "0" "0" "0" "0.1*CONAlpha"  "0"          
## ES    "0" "0" "0" "0" "0" "0.1*ESONAlpha" "0"          
## E     "0" "0" "0" "0" "0" "0"             "0.1*EONBeta"
## I     "0" "0" "0" "0" "0" "0"             "0.1*IONBeta"
## Alpha "0" "0" "0" "0" "0" "0"             "0"          
## Beta  "0" "0" "0" "0" "0" "0"             "0"          
## 
## $S
##       A            C            ES             E            I           
## A     "0.5*AWITHA" "0"          "0"            "0"          "0"         
## C     "0"          "0.5*CWITHC" "0"            "0"          "0"         
## ES    "0"          "0"          "0.5*ESWITHES" "0"          "0"         
## E     "0"          "0"          "0"            "0.5*EWITHE" "0"         
## I     "0"          "0"          "0"            "0"          "0.5*IWITHI"
## Alpha "0"          "0"          "0"            "0"          "0"         
## Beta  "0"          "0"          "0"            "0"          "0"         
##       Alpha             Beta             
## A     "0"               "0"              
## C     "0"               "0"              
## ES    "0"               "0"              
## E     "0"               "0"              
## I     "0"               "0"              
## Alpha "1"               "0*AlphaWITHBeta"
## Beta  "0*AlphaWITHBeta" "1"              
## 
## $F
##    A C ES E I Alpha Beta
## A  1 0  0 0 0     0    0
## C  0 1  0 0 0     0    0
## ES 0 0  1 0 0     0    0
## E  0 0  0 1 0     0    0
## I  0 0  0 0 1     0    0
## 
## $M
##   A C ES E I Alpha Beta
## 1 0 0  0 0 0     0    0

Fixed-effects TSSEM

Stage 1 Analysis

## method="FEM": fixed-effects TSSEM
fixed1 <- tssem1(Digman97.data, Digman97.n, method="FEM")

## summary of the findings
summary(fixed1)
## 
## Call:
## tssem1FEM(Cov = Cov, n = n, cor.analysis = cor.analysis, model.name = model.name, 
##     cluster = cluster, suppressWarnings = suppressWarnings, silent = silent, 
##     run = run)
## 
## Coefficients:
##        Estimate Std.Error z value  Pr(>|z|)    
## S[1,2] 0.363278  0.013368 27.1760 < 2.2e-16 ***
## S[1,3] 0.390375  0.012880 30.3077 < 2.2e-16 ***
## S[1,4] 0.103572  0.015047  6.8830 5.861e-12 ***
## S[1,5] 0.092286  0.015047  6.1331 8.621e-10 ***
## S[2,3] 0.416070  0.012519 33.2345 < 2.2e-16 ***
## S[2,4] 0.135148  0.014776  9.1464 < 2.2e-16 ***
## S[2,5] 0.141431  0.014866  9.5135 < 2.2e-16 ***
## S[3,4] 0.244459  0.014153 17.2724 < 2.2e-16 ***
## S[3,5] 0.138339  0.014834  9.3259 < 2.2e-16 ***
## S[4,5] 0.424566  0.012375 34.3071 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                      Value
## Sample size                      4496.0000
## Chi-square of target model       1505.4443
## DF of target model                130.0000
## p value of target model             0.0000
## Chi-square of independence model 4471.4242
## DF of independence model          140.0000
## RMSEA                               0.1815
## RMSEA lower 95% CI                  0.1736
## RMSEA upper 95% CI                  0.1901
## SRMR                                0.1621
## TLI                                 0.6580
## CFI                                 0.6824
## AIC                              1245.4443
## BIC                               412.0217
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## extract coefficients
coef(fixed1)
##             A         C        ES         E          I
## A  1.00000000 0.3632782 0.3903748 0.1035716 0.09228557
## C  0.36327824 1.0000000 0.4160695 0.1351477 0.14143058
## ES 0.39037483 0.4160695 1.0000000 0.2444593 0.13833895
## E  0.10357155 0.1351477 0.2444593 1.0000000 0.42456626
## I  0.09228557 0.1414306 0.1383390 0.4245663 1.00000000

Stage 2 Analysis

fixed2 <- tssem2(fixed1, RAM=RAM)
summary(fixed2)
## 
## Call:
## wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj), 
##     n = sum(tssem1.obj$n), RAM = RAM, Amatrix = Amatrix, Smatrix = Smatrix, 
##     Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## AONAlpha      0.562800  0.015380 0.532655 0.592944  36.593 < 2.2e-16 ***
## CONAlpha      0.605217  0.015324 0.575183 0.635251  39.495 < 2.2e-16 ***
## EONBeta       0.781413  0.034244 0.714296 0.848529  22.819 < 2.2e-16 ***
## ESONAlpha     0.719201  0.015685 0.688458 0.749943  45.852 < 2.2e-16 ***
## IONBeta       0.551374  0.026031 0.500354 0.602394  21.181 < 2.2e-16 ***
## AlphaWITHBeta 0.362678  0.022384 0.318806 0.406549  16.203 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                4496.0000
## Chi-square of target model                   65.4526
## DF of target model                            4.0000
## p value of target model                       0.0000
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model           3112.7591
## DF of independence model                     10.0000
## RMSEA                                         0.0585
## RMSEA lower 95% CI                            0.0465
## RMSEA upper 95% CI                            0.0713
## SRMR                                          0.0284
## TLI                                           0.9505
## CFI                                           0.9802
## AIC                                          57.4526
## BIC                                          31.8088
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
plot(fixed2, color="green")

Fixed-effects TSSEM with the type of participants as a moderator

Stage 1 analysis

# Display the original study characteristic
table(Digman97.cluster)     
## Digman97.cluster
##   Adolescents      Children Mature adults  Young adults 
##             1             4             6             3
## Younger participants: "Children" and "Adolescents"
## Older participants: "Mature adults"
sample <- ifelse(Digman97.cluster %in% c("Children", "Adolescents"), 
                 yes="Younger participants", no="Older participants")

table(sample)
## sample
##   Older participants Younger participants 
##                    9                    5
## cluster: variable for the analysis with cluster
fixed1.cluster <- tssem1(Digman97.data, Digman97.n, method="FEM", cluster=sample)

summary(fixed1.cluster)
## $`Older participants`
## 
## Call:
## tssem1FEM(Cov = data.cluster[[i]], n = n.cluster[[i]], cor.analysis = cor.analysis, 
##     model.name = model.name, suppressWarnings = suppressWarnings)
## 
## Coefficients:
##        Estimate Std.Error z value  Pr(>|z|)    
## S[1,2] 0.297471  0.015436 19.2716 < 2.2e-16 ***
## S[1,3] 0.370248  0.014532 25.4774 < 2.2e-16 ***
## S[1,4] 0.137694  0.016403  8.3944 < 2.2e-16 ***
## S[1,5] 0.098061  0.016724  5.8637 4.528e-09 ***
## S[2,3] 0.393692  0.014146 27.8306 < 2.2e-16 ***
## S[2,4] 0.183045  0.016055 11.4009 < 2.2e-16 ***
## S[2,5] 0.092774  0.016643  5.5743 2.485e-08 ***
## S[3,4] 0.260753  0.015554 16.7645 < 2.2e-16 ***
## S[3,5] 0.096141  0.016573  5.8009 6.597e-09 ***
## S[4,5] 0.411756  0.013900 29.6224 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                      Value
## Sample size                      3658.0000
## Chi-square of target model        825.9826
## DF of target model                 80.0000
## p value of target model             0.0000
## Chi-square of independence model 3000.9731
## DF of independence model           90.0000
## RMSEA                               0.1515
## RMSEA lower 95% CI                  0.1424
## RMSEA upper 95% CI                  0.1611
## SRMR                                0.1459
## TLI                                 0.7117
## CFI                                 0.7437
## AIC                               665.9826
## BIC                               169.6088
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## 
## $`Younger participants`
## 
## Call:
## tssem1FEM(Cov = data.cluster[[i]], n = n.cluster[[i]], cor.analysis = cor.analysis, 
##     model.name = model.name, suppressWarnings = suppressWarnings)
## 
## Coefficients:
##         Estimate Std.Error z value  Pr(>|z|)    
## S[1,2]  0.604330  0.022125 27.3142 < 2.2e-16 ***
## S[1,3]  0.465536  0.027493 16.9327 < 2.2e-16 ***
## S[1,4] -0.031381  0.035940 -0.8732   0.38258    
## S[1,5]  0.061508  0.034547  1.7804   0.07500 .  
## S[2,3]  0.501432  0.026348 19.0311 < 2.2e-16 ***
## S[2,4] -0.060678  0.034557 -1.7559   0.07911 .  
## S[2,5]  0.320987  0.031064 10.3330 < 2.2e-16 ***
## S[3,4]  0.175437  0.033675  5.2097 1.891e-07 ***
## S[3,5]  0.305149  0.031586  9.6609 < 2.2e-16 ***
## S[4,5]  0.478640  0.026883 17.8045 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                      Value
## Sample size                       838.0000
## Chi-square of target model        346.2810
## DF of target model                 40.0000
## p value of target model             0.0000
## Chi-square of independence model 1470.4511
## DF of independence model           50.0000
## RMSEA                               0.2139
## RMSEA lower 95% CI                  0.1939
## RMSEA upper 95% CI                  0.2355
## SRMR                                0.1411
## TLI                                 0.7305
## CFI                                 0.7844
## AIC                               266.2810
## BIC                                77.0402
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)

Stage 2 analysis

fixed2.cluster <- tssem2(fixed1.cluster, RAM=RAM)

summary(fixed2.cluster)
## $`Older participants`
## 
## Call:
## wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj), 
##     n = sum(tssem1.obj$n), RAM = RAM, Amatrix = Amatrix, Smatrix = Smatrix, 
##     Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## AONAlpha      0.512656  0.018206 0.476973 0.548340  28.158 < 2.2e-16 ***
## CONAlpha      0.549967  0.017945 0.514795 0.585140  30.647 < 2.2e-16 ***
## EONBeta       0.967136  0.058751 0.851986 1.082287  16.462 < 2.2e-16 ***
## ESONAlpha     0.732174  0.018929 0.695074 0.769273  38.680 < 2.2e-16 ***
## IONBeta       0.430649  0.029634 0.372568 0.488730  14.532 < 2.2e-16 ***
## AlphaWITHBeta 0.349236  0.028118 0.294125 0.404346  12.420 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                3658.0000
## Chi-square of target model                   21.9954
## DF of target model                            4.0000
## p value of target model                       0.0002
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model           2273.3179
## DF of independence model                     10.0000
## RMSEA                                         0.0351
## RMSEA lower 95% CI                            0.0217
## RMSEA upper 95% CI                            0.0500
## SRMR                                          0.0160
## TLI                                           0.9801
## CFI                                           0.9920
## AIC                                          13.9954
## BIC                                         -10.8233
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## 
## $`Younger participants`
## 
## Call:
## wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj), 
##     n = sum(tssem1.obj$n), RAM = RAM, Amatrix = Amatrix, Smatrix = Smatrix, 
##     Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##                Estimate Std.Error    lbound    ubound z value Pr(>|z|)    
## AONAlpha       0.747647  0.023880  0.700842  0.794451 31.3081   <2e-16 ***
## CONAlpha       0.911705  0.019864  0.872772  0.950638 45.8969   <2e-16 ***
## EONBeta        0.152563  0.159128 -0.159322  0.464448  0.9587   0.3377    
## ESONAlpha      0.677435  0.025864  0.626743  0.728126 26.1926   <2e-16 ***
## IONBeta        3.283839  3.363262 -3.308033  9.875711  0.9764   0.3289    
## AlphaWITHBeta  0.117257  0.125421 -0.128565  0.363078  0.9349   0.3498    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                 838.0000
## Chi-square of target model                  145.6167
## DF of target model                            4.0000
## p value of target model                       0.0000
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model           2480.2403
## DF of independence model                     10.0000
## RMSEA                                         0.2057
## RMSEA lower 95% CI                            0.1778
## RMSEA upper 95% CI                            0.2350
## SRMR                                          0.1051
## TLI                                           0.8567
## CFI                                           0.9427
## AIC                                         137.6167
## BIC                                         118.6926
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## Setup two plots side-by-side
layout(t(1:2))

## Plot the first group
plot(fixed2.cluster[[1]], col="green")
title("Younger participants")

## Plot the second group
plot(fixed2.cluster[[2]], col="yellow")
title("Older participants")

Random-effects TSSEM

Stage 1 Analysis

## method="REM": Random-effects model
random1 <- tssem1(Digman97.data, Digman97.n, method="REM", RE.type="Diag")

summary(random1)
## 
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues, 
##     "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound, 
##     I2 = I2, model.name = model.name, suppressWarnings = TRUE, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation (robust=FALSE)
## Coefficients:
##                Estimate   Std.Error      lbound      ubound z value  Pr(>|z|)
## Intercept1   0.38971908  0.05429256  0.28330762  0.49613054  7.1781 7.068e-13
## Intercept2   0.43265881  0.04145136  0.35141563  0.51390198 10.4377 < 2.2e-16
## Intercept3   0.04945631  0.06071079 -0.06953466  0.16844728  0.8146   0.41529
## Intercept4   0.09603708  0.04478711  0.00825595  0.18381822  2.1443   0.03201
## Intercept5   0.42724244  0.03911620  0.35057609  0.50390878 10.9224 < 2.2e-16
## Intercept6   0.11929318  0.04106203  0.03881309  0.19977328  2.9052   0.00367
## Intercept7   0.19292424  0.04757962  0.09966990  0.28617858  4.0548 5.018e-05
## Intercept8   0.22690164  0.03240892  0.16338132  0.29042196  7.0012 2.538e-12
## Intercept9   0.18105567  0.04258855  0.09758363  0.26452770  4.2513 2.126e-05
## Intercept10  0.43614968  0.03205960  0.37331402  0.49898535 13.6043 < 2.2e-16
## Tau2_1_1     0.03648372  0.01513055  0.00682839  0.06613905  2.4113   0.01590
## Tau2_2_2     0.01963098  0.00859789  0.00277942  0.03648253  2.2832   0.02242
## Tau2_3_3     0.04571438  0.01952285  0.00745030  0.08397846  2.3416   0.01920
## Tau2_4_4     0.02236122  0.00995083  0.00285794  0.04186449  2.2472   0.02463
## Tau2_5_5     0.01729072  0.00796404  0.00168149  0.03289995  2.1711   0.02992
## Tau2_6_6     0.01815481  0.00895896  0.00059557  0.03571405  2.0264   0.04272
## Tau2_7_7     0.02604881  0.01130266  0.00389602  0.04820161  2.3047   0.02119
## Tau2_8_8     0.00988761  0.00513713 -0.00018098  0.01995619  1.9247   0.05426
## Tau2_9_9     0.01988244  0.00895053  0.00233973  0.03742515  2.2214   0.02633
## Tau2_10_10   0.01049222  0.00501578  0.00066148  0.02032296  2.0918   0.03645
##                
## Intercept1  ***
## Intercept2  ***
## Intercept3     
## Intercept4  *  
## Intercept5  ***
## Intercept6  ** 
## Intercept7  ***
## Intercept8  ***
## Intercept9  ***
## Intercept10 ***
## Tau2_1_1    *  
## Tau2_2_2    *  
## Tau2_3_3    *  
## Tau2_4_4    *  
## Tau2_5_5    *  
## Tau2_6_6    *  
## Tau2_7_7    *  
## Tau2_8_8    .  
## Tau2_9_9    *  
## Tau2_10_10  *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 1220.334
## Degrees of freedom of the Q statistic: 130
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                               Estimate
## Intercept1: I2 (Q statistic)    0.9326
## Intercept2: I2 (Q statistic)    0.8872
## Intercept3: I2 (Q statistic)    0.9325
## Intercept4: I2 (Q statistic)    0.8703
## Intercept5: I2 (Q statistic)    0.8797
## Intercept6: I2 (Q statistic)    0.8480
## Intercept7: I2 (Q statistic)    0.8887
## Intercept8: I2 (Q statistic)    0.7669
## Intercept9: I2 (Q statistic)    0.8590
## Intercept10: I2 (Q statistic)   0.8193
## 
## Number of studies (or clusters): 14
## Number of observed statistics: 140
## Number of estimated parameters: 20
## Degrees of freedom: 120
## -2 log likelihood: -112.4196 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Extract the fixed-effects estimates
(est_fixed <- coef(random1, select="fixed"))
##  Intercept1  Intercept2  Intercept3  Intercept4  Intercept5  Intercept6 
##  0.38971908  0.43265881  0.04945631  0.09603708  0.42724244  0.11929318 
##  Intercept7  Intercept8  Intercept9 Intercept10 
##  0.19292424  0.22690164  0.18105567  0.43614968
## Average correlation matrix
## Convert the estimated vector to a symmetrical matrix
## where the diagonals are fixed at 1 (for a correlation matrix)
averageR <- vec2symMat(est_fixed, diag=FALSE)
dimnames(averageR) <- list(c("A", "C", "ES", "E", "I"),
                           c("A", "C", "ES", "E", "I"))
averageR
##             A         C        ES          E          I
## A  1.00000000 0.3897191 0.4326588 0.04945631 0.09603708
## C  0.38971908 1.0000000 0.4272424 0.11929318 0.19292424
## ES 0.43265881 0.4272424 1.0000000 0.22690164 0.18105567
## E  0.04945631 0.1192932 0.2269016 1.00000000 0.43614968
## I  0.09603708 0.1929242 0.1810557 0.43614968 1.00000000

Stage 2 Analysis

random2 <- tssem2(random1, RAM=RAM)
summary(random2)
## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, RAM = RAM, 
##     Amatrix = Amatrix, Smatrix = Smatrix, Fmatrix = Fmatrix, 
##     diag.constraints = diag.constraints, cor.analysis = cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## AONAlpha      0.569435  0.052425 0.466684 0.672187 10.8619 < 2.2e-16 ***
## CONAlpha      0.590630  0.052649 0.487439 0.693821 11.2182 < 2.2e-16 ***
## EONBeta       0.679955  0.075723 0.531541 0.828370  8.9795 < 2.2e-16 ***
## ESONAlpha     0.760455  0.061963 0.639009 0.881901 12.2726 < 2.2e-16 ***
## IONBeta       0.641842  0.072459 0.499825 0.783859  8.8580 < 2.2e-16 ***
## AlphaWITHBeta 0.377691  0.047402 0.284785 0.470596  7.9679 1.554e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                4496.0000
## Chi-square of target model                    7.8204
## DF of target model                            4.0000
## p value of target model                       0.0984
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model            501.6769
## DF of independence model                     10.0000
## RMSEA                                         0.0146
## RMSEA lower 95% CI                            0.0000
## RMSEA upper 95% CI                            0.0297
## SRMR                                          0.0436
## TLI                                           0.9806
## CFI                                           0.9922
## AIC                                          -0.1796
## BIC                                         -25.8234
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## Plot the parameter estimates
plot(random2, color="green")

An example: A path model

Data preparation

library(readxl)
library(metaSEM)

## Read the Excel file
df <- read_xlsx("Nohe15.xlsx", sheet="Table A1")

head(df)
## # A tibble: 6 × 18
##   Study          N Participants Country Strain RelW1 RelW2 RelS1 RelS2 `% Women`
##   <chr>      <dbl> <chr>        <chr>   <chr>  <dbl> <dbl> <dbl> <dbl>     <dbl>
## 1 Britt & D…   489 Soldiers     United… Depre…  0.94  0.94  0.8   0.8         15
## 2 Demerouti…   335 Employment … the Ne… Exhau…  0.79  0.81  0.85  0.89        70
## 3 Ford (201…   328 Heterogeneo… United… Depre…  0.8   0.81  0.82  0.84        49
## 4 Hammer et…   234 Wives from … United… Depre…  0.91  0.91  0.9   0.9        100
## 5 Hammer et…   234 Husbands fr… United… Depre…  0.9   0.9   0.87  0.87         0
## 6 Innstrand…  2235 Professiona… Norway  Exhau…  0.7   0.7   0.86  0.88        46
## # ℹ 8 more variables: `Publication Status` <chr>, Lag <dbl>, `W1-S2` <dbl>,
## #   `S1-W2` <dbl>, `W1-S1` <dbl>, `W2-S2` <dbl>, `W1-W2` <dbl>, `S1-S2` <dbl>
## Variable names
my.var <- c("W1", "S1", "W2", "S2")

## Split each row as a list
my.list <- split(df, 1:nrow(df))

## Select the correlation coefficients and convert it into a correlation matrix
my.cor <- lapply(my.list, 
                 function(x) {## Convert the correlations into a correlation matrix
                              mat <- vec2symMat(unlist(x[c("W1-S1","W1-W2","W1-S2","S1-W2","S1-S2","W2-S2")]), diag = FALSE)
                              ## Add the dimensions for ease of reference
                              dimnames(mat) <- list(my.var, my.var)
                              ## Return the correlation matrix
                              mat})

## Add the study names for ease of reference
names(my.cor) <- df$Study

head(my.cor)
## $`Britt & Dawson (2005)`
##      W1   S1   W2   S2
## W1 1.00 0.29 0.58 0.22
## S1 0.29 1.00 0.24 0.57
## W2 0.58 0.24 1.00 0.27
## S2 0.22 0.57 0.27 1.00
## 
## $`Demerouti et al. (2004)`
##      W1   S1   W2   S2
## W1 1.00 0.53 0.57 0.41
## S1 0.53 1.00 0.41 0.68
## W2 0.57 0.41 1.00 0.54
## S2 0.41 0.68 0.54 1.00
## 
## $`Ford (2010)`
##      W1   S1   W2   S2
## W1 1.00 0.35 0.75 0.32
## S1 0.35 1.00 0.26 0.74
## W2 0.75 0.26 1.00 0.30
## S2 0.32 0.74 0.30 1.00
## 
## $`Hammer et al. (2005), female subsample`
##      W1   S1   W2   S2
## W1 1.00 0.32 0.57 0.22
## S1 0.32 1.00 0.30 0.43
## W2 0.57 0.30 1.00 0.30
## S2 0.22 0.43 0.30 1.00
## 
## $`Hammer et al. (2005), male subsample`
##      W1   S1   W2   S2
## W1 1.00 0.19 0.54 0.17
## S1 0.19 1.00 0.21 0.60
## W2 0.54 0.21 1.00 0.30
## S2 0.17 0.60 0.30 1.00
## 
## $`Innstrand et al. (2008)`
##      W1   S1   W2   S2
## W1 1.00 0.42 0.63 0.31
## S1 0.42 1.00 0.30 0.62
## W2 0.63 0.30 1.00 0.44
## S2 0.31 0.62 0.44 1.00
## Sample sizes
df$N
##  [1]  489  335  328  234  234 2235   76   94  236  239  239  138  160  151  409
## [16]   78  256  260  600  462  215 1292  470  665  403  153  201  382  130  946
## [31]  730   66
## Display the no. of studies
pattern.na(my.cor, show.na=FALSE)
##    W1 S1 W2 S2
## W1 32 32 32 32
## S1 32 32 32 32
## W2 32 32 32 32
## S2 32 32 32 32
## Display the cumulative sample sizes
pattern.n(my.cor, df$N)
##       W1    S1    W2    S2
## W1 12906 12906 12906 12906
## S1 12906 12906 12906 12906
## W2 12906 12906 12906 12906
## S2 12906 12906 12906 12906

TSSEM

Stage 1 analysis

## Fit the stage one model (fixed-effects model)
fix1 <- tssem1(my.cor, df$N, method="FEM")

summary(fix1)
## 
## Call:
## tssem1FEM(Cov = Cov, n = n, cor.analysis = cor.analysis, model.name = model.name, 
##     cluster = cluster, suppressWarnings = suppressWarnings, silent = silent, 
##     run = run)
## 
## Coefficients:
##         Estimate Std.Error z value  Pr(>|z|)    
## S[1,2] 0.4239505 0.0073295  57.841 < 2.2e-16 ***
## S[1,3] 0.6242536 0.0054290 114.984 < 2.2e-16 ***
## S[1,4] 0.3349734 0.0078888  42.462 < 2.2e-16 ***
## S[2,3] 0.3288615 0.0079183  41.532 < 2.2e-16 ***
## S[2,4] 0.6315353 0.0053669 117.673 < 2.2e-16 ***
## S[3,4] 0.4342405 0.0072612  59.803 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                       Value
## Sample size                      12906.0000
## Chi-square of target model        1524.9269
## DF of target model                 186.0000
## p value of target model              0.0000
## Chi-square of independence model 18047.6811
## DF of independence model           192.0000
## RMSEA                                0.1336
## RMSEA lower 95% CI                   0.1276
## RMSEA upper 95% CI                   0.1400
## SRMR                                 0.1096
## TLI                                  0.9226
## CFI                                  0.9250
## AIC                               1152.9269
## BIC                               -235.6464
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Fit the stage one model (random-effects model)
rand1 <- tssem1(my.cor, df$N, method="REM")

summary(rand1)
## 
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues, 
##     "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound, 
##     I2 = I2, model.name = model.name, suppressWarnings = TRUE, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation (robust=FALSE)
## Coefficients:
##             Estimate Std.Error    lbound    ubound z value  Pr(>|z|)    
## Intercept1 0.3804522 0.0225616 0.3362323 0.4246720 16.8628 < 2.2e-16 ***
## Intercept2 0.6051298 0.0180362 0.5697794 0.6404802 33.5508 < 2.2e-16 ***
## Intercept3 0.3032290 0.0178803 0.2681842 0.3382738 16.9588 < 2.2e-16 ***
## Intercept4 0.3036392 0.0178408 0.2686718 0.3386065 17.0194 < 2.2e-16 ***
## Intercept5 0.6166503 0.0166427 0.5840312 0.6492694 37.0523 < 2.2e-16 ***
## Intercept6 0.3954085 0.0216645 0.3529469 0.4378701 18.2515 < 2.2e-16 ***
## Tau2_1_1   0.0134777 0.0038704 0.0058919 0.0210635  3.4823 0.0004972 ***
## Tau2_2_2   0.0087592 0.0025260 0.0038083 0.0137102  3.4676 0.0005252 ***
## Tau2_3_3   0.0071123 0.0022470 0.0027082 0.0115163  3.1652 0.0015496 ** 
## Tau2_4_4   0.0070585 0.0022121 0.0027229 0.0113941  3.1909 0.0014183 ** 
## Tau2_5_5   0.0072634 0.0021092 0.0031293 0.0113974  3.4436 0.0005740 ***
## Tau2_6_6   0.0122813 0.0034848 0.0054513 0.0191114  3.5243 0.0004246 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 1466.163
## Degrees of freedom of the Q statistic: 186
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.8829
## Intercept2: I2 (Q statistic)   0.8973
## Intercept3: I2 (Q statistic)   0.7743
## Intercept4: I2 (Q statistic)   0.7718
## Intercept5: I2 (Q statistic)   0.8810
## Intercept6: I2 (Q statistic)   0.8748
## 
## Number of studies (or clusters): 32
## Number of observed statistics: 192
## Number of estimated parameters: 12
## Degrees of freedom: 180
## -2 log likelihood: -300.1701 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Extract the fixed-effects estimates
R1 <- coef(rand1, select = "fixed")
R1
## Intercept1 Intercept2 Intercept3 Intercept4 Intercept5 Intercept6 
##  0.3804522  0.6051298  0.3032290  0.3036392  0.6166503  0.3954085
## Convert it to a correlation matrix
R1 <- vec2symMat(R1, diag = FALSE)
R1
##           [,1]      [,2]      [,3]      [,4]
## [1,] 1.0000000 0.3804522 0.6051298 0.3032290
## [2,] 0.3804522 1.0000000 0.3036392 0.6166503
## [3,] 0.6051298 0.3036392 1.0000000 0.3954085
## [4,] 0.3032290 0.6166503 0.3954085 1.0000000
## Add the dimension names for ease of reference
dimnames(R1) <- list(my.var, my.var)
R1
##           W1        S1        W2        S2
## W1 1.0000000 0.3804522 0.6051298 0.3032290
## S1 0.3804522 1.0000000 0.3036392 0.6166503
## W2 0.6051298 0.3036392 1.0000000 0.3954085
## S2 0.3032290 0.6166503 0.3954085 1.0000000

Stage 2 analysis

## Model in lavaan syntax
## The parameter labels are not necessary. But they make the output easier to follow.
model1 <- "W2 ~ w2w*W1 + s2w*S1    # path coefficeients   
           S2 ~ s2s*S1 + w2s*W1
           S1 ~~ w1withs1*W1       # correlation
           S2 ~~ w2withs2*W2
           W1 ~~ 1*W1              # variances of the independent variables are fixed at 1
           S1 ~~ 1*S1"

## Display the model
## See the layout argument in ?semPaths
## layout can be tree, circle, spring, tree2, circle2
plot(model1, color="yellow", layout="spring")

## Convert the above model to RAM specification used by metaSEM.
## We also have to specify how the variables are arranged in the data.
RAM1 <- lavaan2RAM(model1, obs.variables = my.var)
RAM1
## $A
##    W1        S1        W2  S2 
## W1 "0"       "0"       "0" "0"
## S1 "0"       "0"       "0" "0"
## W2 "0.1*w2w" "0.1*s2w" "0" "0"
## S2 "0.1*w2s" "0.1*s2s" "0" "0"
## 
## $S
##    W1           S1           W2             S2            
## W1 "1"          "0*w1withs1" "0"            "0"           
## S1 "0*w1withs1" "1"          "0"            "0"           
## W2 "0"          "0"          "0.5*W2WITHW2" "0*w2withs2"  
## S2 "0"          "0"          "0*w2withs2"   "0.5*S2WITHS2"
## 
## $F
##    W1 S1 W2 S2
## W1  1  0  0  0
## S1  0  1  0  0
## W2  0  0  1  0
## S2  0  0  0  1
## 
## $M
##   W1 S1 W2 S2
## 1  0  0  0  0
## Fit the stage two model
rand2a <- tssem2(rand1, RAM=RAM1)

summary(rand2a)
## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, RAM = RAM, 
##     Amatrix = Amatrix, Smatrix = Smatrix, Fmatrix = Fmatrix, 
##     diag.constraints = diag.constraints, cor.analysis = cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##          Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## s2s      0.586124  0.020790 0.545376 0.626872 28.1926 < 2.2e-16 ***
## w2s      0.080237  0.024842 0.031547 0.128927  3.2299 0.0012385 ** 
## s2w      0.085841  0.024796 0.037242 0.134440  3.4619 0.0005364 ***
## w2w      0.572471  0.022265 0.528834 0.616109 25.7122 < 2.2e-16 ***
## w1withs1 0.380452  0.022562 0.336232 0.424672 16.8628 < 2.2e-16 ***
## w2withs2 0.168885  0.025232 0.119431 0.218338  6.6933 2.182e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                            Value
## Sample size                                12906
## Chi-square of target model                     0
## DF of target model                             0
## p value of target model                        0
## Number of constraints imposed on "Smatrix"     0
## DF manually adjusted                           0
## Chi-square of independence model            3079
## DF of independence model                       6
## RMSEA                                          0
## RMSEA lower 95% CI                             0
## RMSEA upper 95% CI                             0
## SRMR                                           0
## TLI                                         -Inf
## CFI                                            1
## AIC                                            0
## BIC                                            0
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
plot(rand2a, color="green", layout="spring")

Stage 2 analysis with constraints

model2 <- "W2 ~ same*W1 + diff*S1     # regression coefficeients   
           S2 ~ same*S1 + diff*W1
           S1 ~~ w1cs1*W1             # correlation
           S2 ~~ w2cs2*W2
           W1 ~~ 1*W1                 # variances of the independent variables are fixed at 1
           S1 ~~ 1*S1"

## Display the model
plot(model2, color="yellow")

## Convert the above model to RAM specification used by metaSEM.
## We also have to specify how the variables are arranged in the data.
RAM2 <- lavaan2RAM(model2, obs.variables = my.var)

## Fit the stage two model
rand2b <- tssem2(rand1, RAM=RAM2)

summary(rand2b)
## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, RAM = RAM, 
##     Amatrix = Amatrix, Smatrix = Smatrix, Fmatrix = Fmatrix, 
##     diag.constraints = diag.constraints, cor.analysis = cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##       Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## diff  0.082840  0.019712 0.044206 0.121474  4.2026 2.638e-05 ***
## same  0.579844  0.015217 0.550020 0.609669 38.1053 < 2.2e-16 ***
## w1cs1 0.380451  0.022562 0.336231 0.424671 16.8628 < 2.2e-16 ***
## w2cs2 0.168824  0.025241 0.119353 0.218295  6.6886 2.254e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                 Value
## Sample size                                12906.0000
## Chi-square of target model                     0.2247
## DF of target model                             2.0000
## p value of target model                        0.8937
## Number of constraints imposed on "Smatrix"     0.0000
## DF manually adjusted                           0.0000
## Chi-square of independence model            3078.9540
## DF of independence model                       6.0000
## RMSEA                                          0.0000
## RMSEA lower 95% CI                             0.0000
## RMSEA upper 95% CI                             0.0079
## SRMR                                           0.0033
## TLI                                            1.0017
## CFI                                            1.0000
## AIC                                           -3.7753
## BIC                                          -18.7062
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## Comparing nested models with anova()
anova(rand2a, rand2b)
##                 base         comparison ep     minus2LL df       AIC    diffLL
## 1 TSSEM2 Correlation               <NA>  6 1.245611e-19 -6 12.000000        NA
## 2 TSSEM2 Correlation TSSEM2 Correlation  4 2.246694e-01 -4  8.224669 0.2246694
##   diffdf         p
## 1     NA        NA
## 2      2 0.8937451
plot(rand2b, color="green")

OSMASEM: Moderator analysis

## Convert the correlation matrices into a dataframe, which is required in osmasem()
Nohe.df <- Cor2DataFrame(x=my.cor, n=df$N)

## Standardize the moderator "Lag" to improve numerical stability and add it into the dataframe
Nohe.df$data$Lag <- scale(df$Lag)
    
head(Nohe.df$data)
##                                        S1_W1 W2_W1 S2_W1 W2_S1 S2_S1 S2_W2
## Britt...Dawson..2005.                   0.29  0.58  0.22  0.24  0.57  0.27
## Demerouti.et.al...2004.                 0.53  0.57  0.41  0.41  0.68  0.54
## Ford..2010.                             0.35  0.75  0.32  0.26  0.74  0.30
## Hammer.et.al...2005...female.subsample  0.32  0.57  0.22  0.30  0.43  0.30
## Hammer.et.al...2005...male.subsample    0.19  0.54  0.17  0.21  0.60  0.30
## Innstrand.et.al...2008.                 0.42  0.63  0.31  0.30  0.62  0.44
##                                        C(S1_W1 S1_W1) C(W2_W1 S1_W1)
## Britt...Dawson..2005.                     0.001422479   0.0002033635
## Demerouti.et.al...2004.                   0.002076395   0.0002968500
## Ford..2010.                               0.002120708   0.0003031852
## Hammer.et.al...2005...female.subsample    0.002972616   0.0004249776
## Hammer.et.al...2005...male.subsample      0.002972616   0.0004249776
## Innstrand.et.al...2008.                   0.000311227   0.0000444943
##                                        C(S2_W1 S1_W1) C(W2_S1 S1_W1)
## Britt...Dawson..2005.                    0.0008791243   0.0008736611
## Demerouti.et.al...2004.                  0.0012832591   0.0012752845
## Ford..2010.                              0.0013106457   0.0013025009
## Hammer.et.al...2005...female.subsample   0.0018371444   0.0018257278
## Hammer.et.al...2005...male.subsample     0.0018371444   0.0018257278
## Innstrand.et.al...2008.                  0.0001923453   0.0001911500
##                                        C(S2_S1 S1_W1) C(S2_W2 S1_W1)
## Britt...Dawson..2005.                    2.047346e-04   0.0004890894
## Demerouti.et.al...2004.                  2.988513e-04   0.0007139245
## Ford..2010.                              3.052293e-04   0.0007291607
## Hammer.et.al...2005...female.subsample   4.278427e-04   0.0010220714
## Hammer.et.al...2005...male.subsample     4.278427e-04   0.0010220714
## Innstrand.et.al...2008.                  4.479427e-05   0.0001070088
##                                        C(W2_W1 W2_W1) C(S2_W1 W2_W1)
## Britt...Dawson..2005.                    0.0007981106   3.750391e-04
## Demerouti.et.al...2004.                  0.0011650033   5.474451e-04
## Ford..2010.                              0.0011898662   5.591284e-04
## Hammer.et.al...2005...female.subsample   0.0016678466   7.837355e-04
## Hammer.et.al...2005...male.subsample     0.0016678466   7.837355e-04
## Innstrand.et.al...2008.                  0.0001746202   8.205553e-05
##                                        C(W2_S1 W2_W1) C(S2_S1 W2_W1)
## Britt...Dawson..2005.                    3.671018e-04   1.042810e-04
## Demerouti.et.al...2004.                  5.358591e-04   1.522191e-04
## Ford..2010.                              5.472951e-04   1.554677e-04
## Hammer.et.al...2005...female.subsample   7.671487e-04   2.179205e-04
## Hammer.et.al...2005...male.subsample     7.671487e-04   2.179205e-04
## Innstrand.et.al...2008.                  8.031892e-05   2.281583e-05
##                                        C(S2_W2 W2_W1) C(S2_W1 S2_W1)
## Britt...Dawson..2005.                    2.035379e-04    0.001649681
## Demerouti.et.al...2004.                  2.971046e-04    0.002408042
## Ford..2010.                              3.034452e-04    0.002459434
## Hammer.et.al...2005...female.subsample   4.253420e-04    0.003447411
## Hammer.et.al...2005...male.subsample     4.253420e-04    0.003447411
## Innstrand.et.al...2008.                  4.453246e-05    0.000360937
##                                        C(W2_S1 S2_W1) C(S2_S1 S2_W1)
## Britt...Dawson..2005.                    0.0005769098   3.592773e-04
## Demerouti.et.al...2004.                  0.0008421160   5.244375e-04
## Ford..2010.                              0.0008600880   5.356298e-04
## Hammer.et.al...2005...female.subsample   0.0012055935   7.507973e-04
## Hammer.et.al...2005...male.subsample     0.0012055935   7.507973e-04
## Innstrand.et.al...2008.                  0.0001262232   7.860697e-05
##                                        C(S2_W2 S2_W1) C(W2_S1 W2_S1)
## Britt...Dawson..2005.                    0.0008607786   0.0016601406
## Demerouti.et.al...2004.                  0.0012564798   0.0024233097
## Ford..2010.                              0.0012832949   0.0024750266
## Hammer.et.al...2005...female.subsample   0.0017988065   0.0034692681
## Hammer.et.al...2005...male.subsample     0.0017988065   0.0034692681
## Innstrand.et.al...2008.                  0.0001883314   0.0003632254
##                                        C(S2_S1 W2_S1) C(S2_W2 W2_S1)
## Britt...Dawson..2005.                    3.727412e-04   0.0008739022
## Demerouti.et.al...2004.                  5.440909e-04   0.0012756364
## Ford..2010.                              5.557026e-04   0.0013028604
## Hammer.et.al...2005...female.subsample   7.789335e-04   0.0018262316
## Hammer.et.al...2005...male.subsample     7.789335e-04   0.0018262316
## Innstrand.et.al...2008.                  8.155277e-05   0.0001912028
##                                        C(S2_S1 S2_S1) C(S2_W2 S2_S1)
## Britt...Dawson..2005.                    0.0007805638   0.0001951863
## Demerouti.et.al...2004.                  0.0011393902   0.0002849138
## Ford..2010.                              0.0011637064   0.0002909943
## Hammer.et.al...2005...female.subsample   0.0016311782   0.0004078894
## Hammer.et.al...2005...male.subsample     0.0016311782   0.0004078894
## Innstrand.et.al...2008.                  0.0001707811   0.0000427052
##                                        C(S2_W2 S2_W2)        Lag
## Britt...Dawson..2005.                    0.0013981148 -0.6794521
## Demerouti.et.al...2004.                  0.0020408303 -0.7711151
## Ford..2010.                              0.0020843846 -0.8016694
## Hammer.et.al...2005...female.subsample   0.0029217015 -0.1294740
## Hammer.et.al...2005...male.subsample     0.0029217015 -0.1294740
## Innstrand.et.al...2008.                  0.0003058963  0.6038301
## Fit the model without any moderator
## The results are similat to that in TSSEM.
fit0 <- osmasem(model.name="No moderator", RAM=RAM1, data=Nohe.df)
summary(fit0)
## Summary of No moderator 
##  
## free parameters:
##        name  matrix row col    Estimate  Std.Error A    z value     Pr(>|z|)
## 1       w2w      A0  W2  W1  0.57247133 0.02226456    25.712219 0.000000e+00
## 2       w2s      A0  S2  W1  0.08023683 0.02484214     3.229868 1.238476e-03
## 3       s2w      A0  W2  S1  0.08584122 0.02479590     3.461912 5.363533e-04
## 4       s2s      A0  S2  S1  0.58612403 0.02079000    28.192589 0.000000e+00
## 5  w1withs1      S0  S1  W1  0.38045217 0.02256156    16.862847 0.000000e+00
## 6  w2withs2      S0  S2  W2  0.16888459 0.02523192     6.693291 2.182055e-11
## 7    Tau1_1 vecTau1   1   1 -4.30671722 0.28716838   -14.997184 0.000000e+00
## 8    Tau1_2 vecTau1   2   1 -4.73764524 0.28838625   -16.428125 0.000000e+00
## 9    Tau1_3 vecTau1   3   1 -4.94593056 0.31593254   -15.655021 0.000000e+00
## 10   Tau1_4 vecTau1   4   1 -4.95352326 0.31339094   -15.806211 0.000000e+00
## 11   Tau1_5 vecTau1   5   1 -4.92491364 0.29039425   -16.959405 0.000000e+00
## 12   Tau1_6 vecTau1   6   1 -4.39967581 0.28374629   -15.505668 0.000000e+00
## 
## To obtain confidence intervals re-run with intervals=TRUE
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:             12                    180             -300.1701
##    Saturated:             27                    165                    NA
## Independence:             12                    180                    NA
## Number of observations/statistics: 12906/192
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:      -660.1701              -276.1701                -276.1459
## BIC:     -2003.9507              -186.5847                -224.7195
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2024-07-21 15:01:14 
## Wall clock time: 0.0996995 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.21.11 
## Need help?  See help(mxSummary)
plot(fit0)

## Create a matrix to present "Lag" moderator
## We use Lag to moderate the following paths:
## W1 -> W2
## W1 -> S2
## S1 -> W2
## S1 -> S2
A1 <- create.modMatrix(RAM=RAM1, output="A", mod="Lag")
A1
##    W1           S1           W2  S2 
## W1 "0"          "0"          "0" "0"
## S1 "0"          "0"          "0" "0"
## W2 "0*data.Lag" "0*data.Lag" "0" "0"
## S2 "0*data.Lag" "0*data.Lag" "0" "0"
## Fit the model with Lag as a moderator
fit1 <- osmasem(model.name="Lag as a moderator", RAM=RAM1, Ax=A1, data=Nohe.df)
summary(fit1)
## Summary of Lag as a moderator 
##  
## free parameters:
##        name  matrix row col     Estimate  Std.Error A     z value     Pr(>|z|)
## 1       w2w      A0  W2  W1  0.573039774 0.01839766    31.1474347 0.000000e+00
## 2       w2s      A0  S2  W1  0.079844963 0.02419489     3.3000756 9.665879e-04
## 3       s2w      A0  W2  S1  0.085391026 0.02393613     3.5674540 3.604667e-04
## 4       s2s      A0  S2  S1  0.586234249 0.01962561    29.8708837 0.000000e+00
## 5  w1withs1      S0  S1  W1  0.381183701 0.02282277    16.7019024 0.000000e+00
## 6  w2withs2      S0  S2  W2  0.166974817 0.02500439     6.6778197 2.425238e-11
## 7     w2w_1      A1  W2  W1 -0.062015587 0.01850574    -3.3511549 8.047527e-04
## 8     w2s_1      A1  S2  W1 -0.025933718 0.02096724    -1.2368685 2.161359e-01
## 9     s2w_1      A1  W2  S1 -0.002382818 0.02055897    -0.1159016 9.077305e-01
## 10    s2s_1      A1  S2  S1 -0.027809761 0.01974172    -1.4086798 1.589299e-01
## 11   Tau1_1 vecTau1   1   1 -4.276381062 0.28720205   -14.8898001 0.000000e+00
## 12   Tau1_2 vecTau1   2   1 -5.261036068 0.32310515   -16.2827364 0.000000e+00
## 13   Tau1_3 vecTau1   3   1 -5.048387485 0.32015075   -15.7687824 0.000000e+00
## 14   Tau1_4 vecTau1   4   1 -5.039816725 0.31967923   -15.7652303 0.000000e+00
## 15   Tau1_5 vecTau1   5   1 -5.074951110 0.29424551   -17.2473355 0.000000e+00
## 16   Tau1_6 vecTau1   6   1 -4.397726436 0.28288196   -15.5461536 0.000000e+00
## 
## To obtain confidence intervals re-run with intervals=TRUE
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:             16                    176             -323.6921
##    Saturated:             27                    165                    NA
## Independence:             12                    180                    NA
## Number of observations/statistics: 12906/192
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:      -675.6921              -291.6921                -291.6499
## BIC:     -1989.6109              -172.2449                -223.0913
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2024-07-21 15:01:14 
## Wall clock time: 0.2052836 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.21.11 
## Need help?  See help(mxSummary)
## Comparing the models with and without the covariate
anova(fit1, fit0)
##                 base   comparison ep  minus2LL  df       AIC   diffLL diffdf
## 1 Lag as a moderator         <NA> 16 -323.6921 176 -291.6921       NA     NA
## 2 Lag as a moderator No moderator 12 -300.1701 180 -276.1701 23.52199      4
##              p
## 1           NA
## 2 9.957486e-05
## Effect of w2w when Lag is at the mean value
mxEval(w2w, fit1$mx.fit)
##           [,1]
## [1,] 0.5730398
## Effect of w2w when Lag is at +1 SD as Lag is standardized
mxEval(w2w + w2w_1, fit1$mx.fit)
##           [,1]
## [1,] 0.5110242
## Effect of w2w when Lag is at -1 SD as Lag is standardized
mxEval(w2w - w2w_1, fit1$mx.fit)
##           [,1]
## [1,] 0.6350554

An example: A mediation model

Data preparation

library(metaSEM)

## Read the SPSS dataset
my.df <- foreign::read.spss("Schutte21.sav", use.value.labels = TRUE, to.data.frame=TRUE)

## A function to convert rows into a 3x3 correlation matrix
create.matrix <- function(x, type=c(1, 2, 3)) {
  mat <- matrix(NA, ncol=3, nrow=3)
  diag(mat) <- 1
  type <- as.character(type)
  ## Mindfulness, EI, Gratitude
  ## 1: Mindfulness and EI
  ## 2: Mindfulness and Gratitude
  ## 3: EI and Gratitude
  switch(type,
         "1" = mat[1, 2] <- mat[2, 1] <- unlist(x),
         "2" = mat[1, 3] <- mat[3, 1] <- unlist(x),
         "3" = mat[2, 3] <- mat[3, 2] <- unlist(x))
  mat
}

## Convert the data to correlation matrices using the create.matrix().
my.cor <- lapply(split(my.df, seq(nrow(my.df))),
                 function(x, y) create.matrix(x["Effect_size"], x["Type_of_Association"]))

## Variable names
varlist <- c("Mindfulness", "EI", "Gratitude")

## Add the variable names in the dimnames.
my.cor <- lapply(my.cor, function(x) {dimnames(x) <- list(varlist, varlist); x}  )

## Add the study names
names(my.cor) <- my.df$Study_name

## Correlation matrices in the analysis
head(my.cor)
## $`Aras, 2015                                        `
##             Mindfulness   EI Gratitude
## Mindfulness        1.00 0.39        NA
## EI                 0.39 1.00        NA
## Gratitude            NA   NA         1
## 
## $`Bao et al., 2015                                  `
##             Mindfulness   EI Gratitude
## Mindfulness        1.00 0.34        NA
## EI                 0.34 1.00        NA
## Gratitude            NA   NA         1
## 
## $`Bravo et al., 2015                                `
##             Mindfulness   EI Gratitude
## Mindfulness        1.00 0.13        NA
## EI                 0.13 1.00        NA
## Gratitude            NA   NA         1
## 
## $`Brown and Ryan, 2003_Sample A                     `
##             Mindfulness   EI Gratitude
## Mindfulness        1.00 0.46        NA
## EI                 0.46 1.00        NA
## Gratitude            NA   NA         1
## 
## $`Brown and Ryan, 2003_Sample B                     `
##             Mindfulness   EI Gratitude
## Mindfulness        1.00 0.42        NA
## EI                 0.42 1.00        NA
## Gratitude            NA   NA         1
## 
## $`Brown and Ryan, 2003_Sample C                     `
##             Mindfulness   EI Gratitude
## Mindfulness        1.00 0.37        NA
## EI                 0.37 1.00        NA
## Gratitude            NA   NA         1
## Sample sizes
my.n <- my.df$N
my.n
##  [1]  100  380   83  313  327  207 1157  190  256  108  801  365   39   60  250
## [16]  427  470  228  375  234  121  535  156  970  123   99 1392  125  124  200
## [31]  200   72  700  341  313  103  406  299  321  327  200
## Number of studies in each cell
pattern.na(my.cor, show.na = FALSE)
##             Mindfulness EI Gratitude
## Mindfulness          41 26         8
## EI                   26 41         7
## Gratitude             8  7        41
## Total sample sizes in each cell
pattern.n(my.cor, my.n)
##             Mindfulness    EI Gratitude
## Mindfulness       13497  6369      3130
## EI                 6369 13497      3998
## Gratitude          3130  3998     13497

First-stage of analysis

## Stage 1 analysis: find an average correlation matrix
stage1 <- tssem1(my.cor, my.n)
summary(stage1)
## 
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues, 
##     "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound, 
##     I2 = I2, model.name = model.name, suppressWarnings = TRUE, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation (robust=FALSE)
## Coefficients:
##              Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)    
## Intercept1  0.4007062  0.0342472  0.3335830  0.4678294 11.7004 < 2.2e-16 ***
## Intercept2  0.2188135  0.0290679  0.1618414  0.2757855  7.5277 5.174e-14 ***
## Intercept3  0.3123689  0.0416955  0.2306472  0.3940906  7.4917 6.795e-14 ***
## Tau2_1_1    0.0260509  0.0089658  0.0084783  0.0436236  2.9056  0.003666 ** 
## Tau2_2_2    0.0038767  0.0030542 -0.0021094  0.0098628  1.2693  0.204327    
## Tau2_3_3    0.0095986  0.0058884 -0.0019424  0.0211396  1.6301  0.103083    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 262.1055
## Degrees of freedom of the Q statistic: 38
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.9209
## Intercept2: I2 (Q statistic)   0.5784
## Intercept3: I2 (Q statistic)   0.7998
## 
## Number of studies (or clusters): 41
## Number of observed statistics: 41
## Number of estimated parameters: 6
## Degrees of freedom: 35
## -2 log likelihood: -44.57614 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Average correlation matrix
meanR <- vec2symMat(coef(stage1, select = "fixed"), diag = FALSE)
dimnames(meanR) <- list(varlist, varlist)
meanR
##             Mindfulness        EI Gratitude
## Mindfulness   1.0000000 0.4007062 0.2188135
## EI            0.4007062 1.0000000 0.3123689
## Gratitude     0.2188135 0.3123689 1.0000000
## Absolute heterogeneity variance: tau^2
tau2 <- vec2symMat(coef(stage1, select = "random"), diag = FALSE)
dimnames(tau2) <- list(varlist, varlist)
tau2
##             Mindfulness          EI   Gratitude
## Mindfulness  1.00000000 0.026050950 0.003876740
## EI           0.02605095 1.000000000 0.009598596
## Gratitude    0.00387674 0.009598596 1.000000000
## Relative heterogeneity index: I^2
I2 <- vec2symMat(summary(stage1)$I2.values[, "Estimate"], diag = FALSE)
dimnames(I2) <- list(varlist, varlist)
I2
##             Mindfulness       EI Gratitude
## Mindfulness   1.0000000 0.920929 0.5783717
## EI            0.9209290 1.000000 0.7998260
## Gratitude     0.5783717 0.799826 1.0000000

Second-stage of analysis

## Proposed model
model <- "## cp: c prime (c'), the common notation for the direct effect.
          Gratitude ~ cp*Mindfulness + b*EI
          EI ~ a*Mindfulness
          Mindfulness ~~ 1*Mindfulness
          ## Define direct, indirect, and total effects
          Direct := cp
          Indirect := a*b
          Total := a*b + cp"

plot(model, color="yellow")

RAM1 <- lavaan2RAM(model, obs.variables = varlist)
RAM1
## $A
##             Mindfulness EI      Gratitude
## Mindfulness "0"         "0"     "0"      
## EI          "0.1*a"     "0"     "0"      
## Gratitude   "0.1*cp"    "0.1*b" "0"      
## 
## $S
##             Mindfulness EI             Gratitude                   
## Mindfulness "1"         "0"            "0"                         
## EI          "0"         "0.5*EIWITHEI" "0"                         
## Gratitude   "0"         "0"            "0.5*GratitudeWITHGratitude"
## 
## $F
##             Mindfulness EI Gratitude
## Mindfulness           1  0         0
## EI                    0  1         0
## Gratitude             0  0         1
## 
## $M
##   Mindfulness EI Gratitude
## 1           0  0         0
## 
## $mxalgebras
## $mxalgebras$Direct
## mxAlgebra 'Direct' 
## $formula:  cp 
## $result: (not yet computed) <0 x 0 matrix>
## dimnames: NULL
## 
## $mxalgebras$Indirect
## mxAlgebra 'Indirect' 
## $formula:  a * b 
## $result: (not yet computed) <0 x 0 matrix>
## dimnames: NULL
## 
## $mxalgebras$Total
## mxAlgebra 'Total' 
## $formula:  a * b + cp 
## $result: (not yet computed) <0 x 0 matrix>
## dimnames: NULL
## Stage 2 analysis: fit the path model
## Likelihood-based CI: intervals.type = "LB"
stage2 <- tssem2(stage1, RAM=RAM1, intervals.type = "LB")
summary(stage2)
## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, RAM = RAM, 
##     Amatrix = Amatrix, Smatrix = Smatrix, Fmatrix = Fmatrix, 
##     diag.constraints = diag.constraints, cor.analysis = cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: Likelihood-based statistic
## Coefficients:
##    Estimate Std.Error   lbound   ubound z value Pr(>|z|)
## a  0.400706        NA 0.333466 0.467842      NA       NA
## b  0.267667        NA 0.166085 0.369322      NA       NA
## cp 0.111558        NA 0.028737 0.190505      NA       NA
## 
## mxAlgebras objects (and their 95% likelihood-based CIs):
##                   lbound  Estimate    ubound
## Direct[1,1]   0.02873693 0.1115576 0.1905046
## Indirect[1,1] 0.06580411 0.1072559 0.1559166
## Total[1,1]    0.16178591 0.2188135 0.2760576
## 
## Goodness-of-fit indices:
##                                               Value
## Sample size                                13497.00
## Chi-square of target model                     0.00
## DF of target model                             0.00
## p value of target model                        0.00
## Number of constraints imposed on "Smatrix"     0.00
## DF manually adjusted                           0.00
## Chi-square of independence model             249.69
## DF of independence model                       3.00
## RMSEA                                          0.00
## RMSEA lower 95% CI                             0.00
## RMSEA upper 95% CI                             0.00
## SRMR                                           0.00
## TLI                                            -Inf
## CFI                                            1.00
## AIC                                            0.00
## BIC                                            0.00
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
plot(stage2, color="green")

Testing the direct effect equals the indirect effect (advanced topic!!!)

## Proposed model
model <- "Gratitude ~ cp*Mindfulness + b*EI
          EI ~ a*Mindfulness
          Mindfulness ~~ 1*Mindfulness
          ## Define direct, indirect, and total effects
          Direct := cp
          Indirect := a*b
          Total := a*b + cp
          ## Apply a constraint on the direct and indirect effects
          cp == a*b"

RAM2 <- lavaan2RAM(model, obs.variables = varlist)
RAM2
## $A
##             Mindfulness EI      Gratitude
## Mindfulness "0"         "0"     "0"      
## EI          "0.1*a"     "0"     "0"      
## Gratitude   "0.1*cp"    "0.1*b" "0"      
## 
## $S
##             Mindfulness EI             Gratitude                   
## Mindfulness "1"         "0"            "0"                         
## EI          "0"         "0.5*EIWITHEI" "0"                         
## Gratitude   "0"         "0"            "0.5*GratitudeWITHGratitude"
## 
## $F
##             Mindfulness EI Gratitude
## Mindfulness           1  0         0
## EI                    0  1         0
## Gratitude             0  0         1
## 
## $M
##   Mindfulness EI Gratitude
## 1           0  0         0
## 
## $mxalgebras
## $mxalgebras$Direct
## mxAlgebra 'Direct' 
## $formula:  cp 
## $result: (not yet computed) <0 x 0 matrix>
## dimnames: NULL
## 
## $mxalgebras$Indirect
## mxAlgebra 'Indirect' 
## $formula:  a * b 
## $result: (not yet computed) <0 x 0 matrix>
## dimnames: NULL
## 
## $mxalgebras$Total
## mxAlgebra 'Total' 
## $formula:  a * b + cp 
## $result: (not yet computed) <0 x 0 matrix>
## dimnames: NULL
## 
## $mxalgebras$constraint1
## MxConstraint 'constraint1' 
## $formula:  cp == a * b
stage2b <- tssem2(stage1, RAM=RAM2, intervals.type = "LB")
summary(stage2b)
## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, RAM = RAM, 
##     Amatrix = Amatrix, Smatrix = Smatrix, Fmatrix = Fmatrix, 
##     diag.constraints = diag.constraints, cor.analysis = cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: Likelihood-based statistic
## Coefficients:
##    Estimate Std.Error   lbound   ubound z value Pr(>|z|)
## a  0.401476        NA 0.338588 0.465932      NA       NA
## b  0.270746        NA 0.215989 0.330137      NA       NA
## cp 0.108698        NA 0.087826 0.130046      NA       NA
## 
## mxAlgebras objects (and their 95% likelihood-based CIs):
##                   lbound Estimate    ubound
## Direct[1,1]   0.08782637 0.108698 0.1300457
## Indirect[1,1] 0.08782431 0.108698 0.1299560
## Total[1,1]    0.17576707 0.217396 0.2599277
## 
## Goodness-of-fit indices:
##                                                 Value
## Sample size                                13497.0000
## Chi-square of target model                     0.0052
## DF of target model                             0.0000
## p value of target model                        0.0000
## Number of constraints imposed on "Smatrix"     0.0000
## DF manually adjusted                           0.0000
## Chi-square of independence model             249.6905
## DF of independence model                       3.0000
## RMSEA                                          0.0000
## RMSEA lower 95% CI                             0.0000
## RMSEA upper 95% CI                             0.0000
## SRMR                                           0.0015
## TLI                                              -Inf
## CFI                                            1.0000
## AIC                                            0.0052
## BIC                                            0.0052
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## Compare the models with and without c=a*b
anova(stage2, stage2b)
##                 base         comparison ep     minus2LL df      AIC      diffLL
## 1 TSSEM2 Correlation               <NA>  3 2.893374e-24 -3 6.000000          NA
## 2 TSSEM2 Correlation TSSEM2 Correlation  3 5.222545e-03 -2 6.005223 0.005222545
##   diffdf         p
## 1     NA        NA
## 2      1 0.9423893

Tips in MASEM

It takes a long time to run large models!

library(metaSEM)

## Use all cores-2, i.e., keep 2 cores for other things.
mxOption(key='Number of Threads', value=parallel::detectCores()-2)

How do I know if my correlation matrices are “proper”?

## X1 is okay
X1 <- matrix(c(1, 0.8, 0.8, 1), nrow=2)
X1
##      [,1] [,2]
## [1,]  1.0  0.8
## [2,]  0.8  1.0
is.pd(X1)
## [1] TRUE
## X2 is not okay as the correlation is 1.2.
X2 <- matrix(c(1, 1.2, 1.2, 1), nrow=2)
X2
##      [,1] [,2]
## [1,]  1.0  1.2
## [2,]  1.2  1.0
is.pd(X2)
## [1] FALSE
## X3 is not okay as the correlation matrix is not symmetric.
X3 <- matrix(c(1, 0.6, 0.8, 1), nrow=2)
X3
##      [,1] [,2]
## [1,]  1.0  0.8
## [2,]  0.6  1.0
isSymmetric(X3)
## [1] FALSE

What should I do when the solution is non-convergent?

OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values indicate problems.)
OpenMx status1: 5 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)
Warning message:
In print.summary.meta(x) :
  OpenMx status1 is neither 0 or 1. You are advised to 'rerun' it again.

## Suppose the results are stored in an object called "fit"
fit <- rerun(fit)
summary(fit)
pattern.na(Hunter83$data, show.na=FALSE)
##               Ability Job_knowledge Work_sample Supervisor
## Ability            13            11          11         12
## Job_knowledge      11            12          10         11
## Work_sample        11            10          12         11
## Supervisor         12            11          11         13

Exercises

Further resources

Thank you!

https://m.xkcd.com/1447/

  1. Browne, M. W., & Cudeck, R. (1993). Alternative ways of assessing model fit. In K. A. Bollen & J. S. Long (Eds.), Testing Structural Equation Models (pp. 136-162). Newbury Park, CA: Sage.↩︎

  2. Hu, L., & Bentler, P. M. (1999). Cutoff criteria for fit indices in covariance structure analysis: Conventional criteria versus new alternatives. Structural Equation Modeling, 6, 1-55.↩︎

  3. Marsh, H. W., & Hau, K. T., & Wen, Z. (2004). In search of golden rules: Comment on hypothesis-testing approaches to setting cutoff values for fit indices and dangers in overgeneralizing Hu and Bentler’s (1999) findings. Structural Equation Modeling, 11, 320-342.↩︎

  4. Guyatt, G. H., Sackett, D. L., Sinclair, J. C., Hayward, R., Cook, D. J., Cook, R. J., Bass, E., Gerstein, H., Haynes, B., Holbrook, A., Jaeschke, R., Laupacls, A., Moyer, V., & Wilson, M. (1995). Users’ Guides to the Medical Literature: IX. A Method for Grading Health Care Recommendations. JAMA, 274(22), 1800–1804. https://doi.org/10.1001/jama.1995.03530220066035↩︎

  5. Schmidt, F. L., & Hunter, J. E. (2015). Methods of meta-analysis: Correcting error and bias in research findings (3rd ed.). Sage.↩︎

  6. Hedges, L. V., & Olkin, I. (1985). Statistical methods for meta-analysis. Orlando, FL: Academic Press.↩︎

  7. Hedges, L. V., & Vevea, J. L. (1998). Fixed- and random-effects models in meta-analysis. Psychological Methods, 3(4), 486–504. https://doi.org/10.1037/1082-989X.3.4.486↩︎

  8. Hedges, L. V., & Vevea, J. L. (1998). Fixed- and random-effects models in meta-analysis. Psychological Methods, 3, 486-504.↩︎

  9. Borenstein, M., Hedges, L. V., Higgins, J. P. T., & Rothstein, H. R. (2009). Introduction to meta-analysis. John Wiley & Sons.↩︎

  10. Higgins, J. P. T., & Thompson, S. G. (2002). Quantifying heterogeneity in a meta-analysis. Statistics in Medicine, 21, 1539-1558.↩︎

  11. Borenstein, M., Higgins, J. P. T., Hedges, L. V., & Rothstein, H. R. (2017). Basics of meta-analysis: I2 is not an absolute measure of heterogeneity. Research Synthesis Methods, 8(1), 5–18. https://doi.org/10.1002/jrsm.1230↩︎

  12. Rubio-Aparicio, M., López-López, J. A., Viechtbauer, W., Marín-Martínez, F., Botella, J., & Sánchez-Meca, J. (2020). Testing categorical moderators in mixed-effects meta-analysis in the presence of heteroscedasticity. The Journal of Experimental Education, 88(2), 288–310. https://doi.org/10.1080/00220973.2018.1561404↩︎

  13. Konstantopoulos, S., & Hedges, L. V. (2004). Meta-analysis. In D. Kaplan (Ed.), Handbook of quantitative methodology for the social sciences (pp. 281-297). New York: Sage.↩︎

  14. Cheung, M. W.-L. (2021). Meta-analytic structural equation modeling. In Oxford Research Encyclopedia of Business and Management. Oxford University Press. https://doi.org/10.1093/acrefore/9780190224851.013.225↩︎

  15. National Research Council (1992). Combining information: Statistical issues and opportunities for research. Washington, D.C.: National Academy Press.↩︎

  16. Brown, S. P., & Stayman, D. M. (1992). Antecedents and consequences of attitude toward the ad: A meta-analysis. Journal of Consumer Research, 19, 34-51.↩︎

  17. Premack, S. L., & Hunter, J. E. (1988). Individual unionization decisions. Psychological Bulletin, 103, 223-234.↩︎

  18. Norton, S., Cosco, T., Doyle, F., Done, J., & Sacker, A. (2013). The Hospital Anxiety and Depression Scale: A meta confirmatory factor analysis. Journal of Psychosomatic Research, 74(1), 74-81.↩︎

  19. Viswesvaran, C., & Ones, D. S. (1995). Theory testing: Combining psychometric meta-analysis and structural equations modeling. Personnel Psychology, 48(4), 865-885.↩︎

  20. Becker, B. J. (1992). Using results from replicated studies to estimate linear models. Journal of Educational Statistics, 17(4), 341-362.↩︎

  21. Cheung, M. W.-L. (2014). Fixed- and random-effects meta-analytic structural equation modeling: Examples and analyses in R. Behavior Research Methods, 46(1), 29-40. http://doi.org/10.3758/s13428-013-0361-y↩︎

  22. Cheung, M. W.-L., & Chan, W. (2005). Meta-analytic structural equation modeling: A two-stage approach. Psychological Methods, 10(1), 40-64. http://doi.org/10.1037/1082-989X.10.1.40↩︎

  23. Jak, S., & Cheung, M. W.-L. (2020). Meta-analytic structural equation modeling with moderating effects on SEM parameters. Psychological Methods, 25(4), 430–455. https://doi.org/10.1037/met0000245↩︎

  24. Cudeck, R. (1989). Analysis of correlation matrices using covariance structure models. Psychological Bulletin, 105(2), 317-327.↩︎

  25. Cudeck, R. (1989). Analysis of correlation matrices using covariance structure models. Psychological Bulletin, 105(2), 317–327. https://doi.org/10.1037/0033-2909.105.2.317↩︎

  26. Jak, S., & Cheung, M. W.-L. (2020). Meta-analytic structural equation modeling with moderating effects on SEM parameters. Psychological Methods, 25(4), 430–455. https://doi.org/10.1037/met0000245↩︎

  27. Jak, S., & Cheung, M. W. (2022, March 16). Can findings from meta-analytic structural equation modeling in management and organizational psychology be trusted?. https://psyarxiv.com/b3qvn.↩︎

  28. Cheung, M. W.-L., & Chan, W. (2005). Meta-analytic structural equation modeling: A two-stage approach. Psychological Methods, 10(1), 40–64. https://doi.org/10.1037/1082-989X.10.1.40↩︎

  29. Cheung, M. W.-L. (2014). Fixed- and random-effects meta-analytic structural equation modeling: Examples and analyses in R. Behavior Research Methods, 46(1), 29–40. https://doi.org/10.3758/s13428-013-0361-y↩︎

  30. Digman, J. M. (1997). Higher-order factors of the Big Five. Journal of Personality and Social Psychology, 73(6), 1246-1256. http://doi.org/10.1037/0022-3514.73.6.1246↩︎

  31. Nohe, C., Meier, L. L., Sonntag, K., & Michel, A. (2015). The chicken or the egg? A meta-analysis of panel studies of the relationship between work–family conflict and strain. Journal of Applied Psychology, 100(2), 522–536. https://doi.org/10.1037/a0038012↩︎

  32. Schutte, N. S., Keng, S.-L., & Cheung, M. W.-L. (2021). Emotional intelligence (EI) mediates the connection between mindfulness and gratitude: A meta-analytic structural equation modeling study. Mindfulness, 12(11), 2613–2623. https://doi.org/10.1007/s12671-021-01725-2↩︎

  33. Cheung, M. W.-L. (2021). Meta-analytic structural equation modeling. In Oxford Research Encyclopedia of Business and Management. Oxford University Press. https://doi.org/10.1093/acrefore/9780190224851.013.225↩︎

  34. Cheung, M. W.-L. (2022). Synthesizing indirect effects in mediation models with meta-analytic methods. Alcohol and Alcoholism, 57(1), 5–15. https://doi.org/10.1093/alcalc/agab044↩︎

  35. Cheung, M. W.-L., & Hafdahl, A. R. (2016). Special issue on meta-analytic structural equation modeling: Introduction from the guest editors. Research Synthesis Methods, 7(2), 112–120. https://doi.org/10.1002/jrsm.1212↩︎

  36. Cheung, M. W.-L., & Hong, R. Y. (2017). Applications of meta-analytic structural equation modelling in health psychology: Examples, issues, and recommendations. Health Psychology Review, 11(3), 265–279. https://doi.org/10.1080/17437199.2017.1343678↩︎

  37. Valentine, J. C., Cheung, M. W.-L., Smith, E. J., Alexander, O., Hatton, J. M., Hong, R. Y., Huckaby, L. T., Patton, S. C., Pössel, P., & Seely, H. D. (2022). A primer on meta-analytic structural equation modeling: The case of depression. Prevention Science, 23(3), 346–365. https://doi.org/10.1007/s11121-021-01298-5↩︎